Synchronizing distant nodes: a universal classification of networks 
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Stability of synchronization in delay-coupled networks of identical units generally depends in a 
complicated way on the coupling topology. We show that for large coupling delays synchronizability 
relates in a simple way to the spectral properties of the network topology. The master stability 
function used to determine stability of synchronous solutions has a universal structure in the limit 
of large delay: it is rotationally symmetric around the origin and increases monotonically with the 
radius in the complex plane. This allows a universal classification of networks with respect to their 
synchronization properties and solves the problem of complete synchronization in networks with 
strongly delayed coupling. 
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Synchronization phenomena in networks are of great 
importance [T] in many areas. Chaos synchronization of 
lasers, for instance, may lead to new secure communi- 
cation schemes |2_. The synchronization of neurons is 
believed to play a crucial role in the brain under nor- 
mal conditions, for instance in the context of cognition 
and learning [3 , and under pathological conditions such 
as Parkinson's disease [I]. Time delay effects are a key 
issue in realistic networks. For example, the finite propa- 
gation time of light between coupled semiconductor lasers 
[5 significantly influences the dynamics. Similar effects 
occur in neuronal [6. and biological [7] networks. 

To determine the stability of a synchronized state in a 
network of identical units, a powerful method has been 
developed j8j, i.e., the master stability function (MSF). 
Recent works (HIEI] have started to investigate the MSF 
for networks with coupling delays and found that the 
MSF depends non-trivially on delay times. 

In this work we show that in the limit of large cou- 
pling delays the MSF has a very simple structure. This 
solves the problem of complete zero-lag synchronization 
for networks with large coupling delay. After briefly in- 
troducing the notion of the MSF, we demonstrate the 
implications for large coupling delays based on a scaling 
theory [1 lj . This allows us to describe the synchronizabil- 
ity of networks with strongly delayed coupling depending 
on the type of node dynamics and spectral properties of 
the network topology. For example, as recently conjec- 
tured [10 , networks for which the trajectory of an uncou- 
pled unit is also a solution of the network cannot exhibit 
chaos synchronization for large coupling delay. The re- 
sults presented here confirm and generalize these previous 
findings. 

Consider a system of N identical units connected in 
a network with a coupling delay r [HJ] (x l S M. d , i — 
1,...,N) 

±<(i) = / + Y^ =1 9ijh [x*(t - t)] . (1) 

Here, gij is the real-valued coupling matrix, which deter- 
mines the topology and the strength of each link in the 
network, / is a (non-linear) function describing the dy- 
namics of an isolated unit, and ft, is a possibly non-linear 



coupling function. To allow for an invariant synchroniza- 
tion manifold (SM), the row sum a = Ylj=x 9ij °f the 
matrix has to be the same for each row i [5]. The sta- 
bility of the synchronized solution is then governed by 
the MSF and the eigenvalues of the coupling matrix g^ . 
The MSF is defined as the maximum Lyapunov exponent 
Ama^re 4 ^) as a function of the complex argument re*^ 
arising from the variational equation 

i{t) = Df[x(t)} £(t) + re**Dh[x(t - r)] £(t - r), 

where x(t) is given by the dynamics within the SM. The 
synchronized state is stable for a given coupling topology 
if the MSF is negative at all transversal eigenvalues 7^ 
of the coupling matrix (A max (7fe) < 0). Here, transversal 
eigenvalue refers to all eigenvalues except for the eigen- 
value a associated to perturbations within the SM with 
corresponding eigenvector (1, 1, 1). 

We will now restrict our analysis to maps [TU], but all 
ingredients of our argument are also valid for flows. For 
delay-coupled maps the dynamics in the SM is governed 
by the equation Xk+i = f(%k) + &h(xk-r) with r G N 
and Xk £ C d or £ M. d and the MSF is calculated for fixed 
a from 

£ fe +i = A k Z k + re^B k ^ T (2) 

with matrices — Df(xk) and = Dh(xk- T )- 

Note that when the delay is changed the dynamics in 
the SM changes, too. Hence, we are not able to make pre- 
dictions about what happens as t is changed. However, 
at a fixed large value of the delay time r we can analyze 
the Lyapunov exponents arising from different values of 
re 1 ^ in Eq. We do this in the following steps: first we 
analyze the two simpler cases when the dynamics in the 
SM is a fixed point (FP) or a periodic orbit (PO). Then, 
to expand the results to chaotic dynamics in the SM, we 
use the fact that POs are dense in a chaotic attractor. 

For FPs and POs of delay differential equations a scal- 
ing theory for the eigenvalues or Floquet exponents in 
the limit of large delay [TT] shows that the spectrum con- 
sists in both cases of two parts: a strongly unstable part 
arising from unstable eigenvalues of the system without 
delay and a pseudo-continuous spectrum for which the 



real parts of the eigenvalues approach zero in the limit of 
large delay. This scaling theory has been developed for 
flows; to prove our statements we will extend this theory 
to maps. 

Fixed point - Let us first consider the case of a FP in 
the SM, for which A — Ak and B = Bk are constant. 
Making the ansatz = z fe £oi we find an equation for the 
multipliers z 



detL4- zl + re^B ; 



0. 



(3) 



where I denotes the identity matrix. 

For the strongly unstable spectrum we suppose there 
is a solution with \z\ > 1. Then in the limit of t — > oo 
Eq. ^ becomes detL4 - zl] = 0. Thus in the limit of 
large delay the eigenvalues z of A with \z\ > 1 are also 
solutions of Eq. ^ and vice versa. 

We are now interested in the pseudo-continuous spec- 
trum, i. e., the solutions with \z\ w 1 in the limit of large 
t. We make the ansatz z = (1 + S/r)e lul . In the limit 
t — y oo we have (1 + 5/t)~' 
and Eq. ^ becomes 



e- 5 and (1 + S/t) ~ 1, 



det[A - Ie* 



>B] =0 



(4) 



with (f> — lot. As we will show below, uj as well as the 
parameter <p take on any (arbitrarily dense) values in 
[— 7T, 7r]. From this it is clear that the phase ip in the 
variational equation does not change 5, i.e., the MSF is 
invariant under phase shifts (rotations) and its value only 
depends on r. 

Equation Q is a polynomial in /i = re~ s e l ^~^ for 
which the roots can be calculated. For example, if B is 
invertible, the roots fi are the eigenvalues of the matrix 
— B^ 1 {A — Ie lu ). In general, each root n is a function of 
lj and one can find the branches 5(lo) = — In |/i(w)| + In r 
from the definition of /i. The function /i(w) can admit 
the zero value at some point ljq, i.e., //(wo) — 0, in the 
case when the matrix A has an eigenvalue with \z\ = 1. 
Indeed, as follows from Eq. Q, for fi = 0, U) = ojq and 
det B ^ we have det [A - Ie iuJ °] = det[A - Iz] = 0. In 
all other cases, with det and \z\ ^ 1, the function 

|/Lt(w)| is bounded < /io < 1/^(^)1 < Ml- 

If there are no strongly unstable eigenvalues, the sign 
of 5 determines the stability in the limit of large r, since 
|z|~|1 + <5/t|. It is clear that 6 increases monotonically 
with increasing r and in particular 6 is negative for small 
r and positive for large r. Thus there is a critical radius 
ro for which the first eigenvalue branch becomes unstable 
(6 > 0) and thus the MSF changes sign. 

Note that we have obtained the function 5(uj) on which 
the solutions lie in the limit of large r but not yet the 
exact values of lj. These values can be calculated from 
the expression /i(w) = re _<5 ^'e l ^~ WT - ) , which implies 



Arg n(u) = ip — lot + 2nk 



(5) 



for any integer k. Since /i(w) is a known root of Eq. Q, 
Eq. (5) can be considered as a transcendental equation for 
determining the solutions lu = uj^. In particular, Eq. ^ 
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Figure 1: (Color online) (a) Pseudo-continuous spectrum 
S(uj) (lines) and location of the exact roots (crosses) for a 
one-dimensional complex map for r = 3.3 > ro = 3 and 
r = 2.7 < r = 3. Parameters: A = 0.4, B = 0.2, tp = 0, 
r = 30. (b) Contour line A max = of the MSF for coupled 
semiconductor lasers according to Eq. ( 12 1 for delay times 
t = 1 (solid), 8 (dashed), 20 (dash-dotted)7and 1000 (dotted). 



implies that the distance between neighboring solutions 
uj k and w fe _i 

u k - ojk^x = [Arg^(w fc _i) - Arg^(w fe )]/r + 2tt/t 
= 2it/t + 0(1/t 2 ) 

is proportional to 1/r and the curve S(ui) is filled densely 
with equally spaced roots as r —¥ oo. 

For illustration, consider the simple case of a one- 
dimensional complex map with A, B G C with \A\ < 
1. In this case we can explicitly calculate the pseudo- 
continuous spectrum 6(w) — \n(\rB\/\A — e IW |), which 
is depicted in Fig. [lji. For r < (1 - |A|)/|S| all the 
eigenvalues approach \z\ = 1 from the stable side and 
for r > (1 — \A\)/\B\ there are always weakly unsta- 
ble eigenvalues. Thus the critical radius is given by 
r = (l-|A|)/|B|. 

Periodic orbit - Now consider the variational Eq. ^ 
with Ak and Bk being periodic in k with period T, cor- 
responding to a PO in the SM. We consider the case of 
large delay, i.e., r T. Making a Floquet-like ansatz 
£k = z k qk, where qk is T-periodic, we find 



zq k +i = A k q k + re^Bk z" 



(6) 



with n = rmodT e {0, 1, . . . , T - 1}. 

For the strongly unstable spectrum again suppose 
there is a solution with \z\ > 1, then in the limit r — > oo 
the term z~ T vanishes and we find 

zq k+ i = A k q k . (7) 

Using the periodicity of q kl Eq. ^ implies det[z T — 
rifc=i Ak] — 0, where z T is a Floquet multiplier of the 
system £, k +i — Ak£,k without delay. Hence, if z T is a Flo- 
quet multiplier of Eq. ([7]), with \z\ > 1, then in the limit 
t — > oo it is also a solution of Eq. Q and vice versa. 

For the pseudo-continuous spectrum we again make 
the ansatz z = ( \ + 8 / T)e luJ . Taking the limit r — ► oo Eq. 



^ becomes 

e w q k+1 = A k q k + re- s e^-^B h q k ^ n (8) 

with <f> — u>t. Thus one has to solve 

[e iu J + A + nB}q = 0, (9) 

where fi — re" s e 1 ^~^ and q = (qi, . . . , qj-). The ma- 
trices J, A, and B follow from Eq. taking into 
account the periodicity of A k , B k , and q k , e.g., A — 
diagj^!, . . . , Ax}- Taking the determinant of the en- 
tire matrix in Eq. ^ results in a polynomial in /i = 
re~ s e l ^~^ (of maximum order Td). Again, the roots 
/I are functions of u> and we can calculate the branches 
5(u>) = — In + lnr, where ip and <f> drop out. As in 
the case of FPs, one can show that the function |/z(w)| is 
bounded < < < Mi unless the instantaneous 
system has a Floquet multiplier z with \z\ = 1. Note 
that for the FP case as well as for the PO case, one can 
show that the discussed strongly unstable and pseudo- 
continuous spectrum constitute the entire spectrum. 

We have found the same structure of the MSF for a 
PO in the SM: The MSF is rotationally symmetric about 
the origin in the complex plane. If without feedback (r — 
0) the MSF is positive, then it is a positive constant in 
the limit of large delay. Otherwise it is a monotonically 
increasing function of r and it changes sign at a critical 
radius ro- 

Chaotic dynamics - Every chaotic attractor embeds 
an infinite number of unstable periodic orbits (UPOs). 
It is well known that the characteristic properties of the 
chaotic system can be described in terms of these UPOs. 
One of the most important examples is the natural mea- 
sure of the chaotic attractor which is concentrated at the 
UPOs and can in fact be expressed in terms of the orbit's 
Floquet multipliers [T2|[T3]. 

Lyapunov exponents arising from variational equations 
such as Eq. Q have been discussed in the framework of 
PO theory [14 , too. In particular it has been shown fT5] 
that a chaotic attractor in an invariant manifold loses 
its transversal stability in a blow-out bifurcation when 
the transversely unstable orbits outweigh the transversely 
stable orbits. To be precise, we divide the orbits into 
these two groups and define [T5| the transversely stable 
weight Afn and the unstable weight Aj, as 



Table I: Stability of chaotic and non-chaotic synchronized 
solutions for the three types of networks 



' rf 7 — 



(10) 



where the sum goes over all Nj. transversely unstable and 
TVj. transversely stable orbits with period T (or factors 
of T), respectively. Here, ht{j) is the weight of the jth 
orbit, corresponding to the natural measure of a typical 
trajectory in the neighborhood of the jth orbit and Xx(j) 
is the transversal Lyapunov exponent of this jth orbit. 
The weight of a PO is inversely proportional to the prod- 
uct of its unstable Floquet multipliers [12) . The attractor 
is transversely unstable if and only if in the limit of large 
T 





chaotic dynamics 
in the SM (r < \a\) 


PO or FP in the SM 

(kl < ro) 


(a) 


7max < |c| 


synchr. stable if 

7max < To 


synchr. stable 


(b) 


7max = |c| 


synchr. unstable 


synchr. stable 


(c) 


7max > O" 


synchr. unstable 


synchr. stable 

if 7max < To 



We now draw the connection to the scaling theory for 
large r. Starting from r — (no feedback), the transver- 
sal Lyapunov exponents Xt(j) of each orbit can only in- 
crease with increasing r, as shown above, and the weights 
[ir (J) are not changed. In particular for large enough 
r any orbit becomes transversely unstable: either it is 
already unstable for r = and thus remains unstable, 
or the pseudo-continuous spectrum goes to zero and for 
large r it does so from the unstable side. Thus there ex- 



ists a minimum radius ro for which the condition (111 
on the weights is fulfilled. Note that since we consider 
the limit r — > oo we can evaluate Eq. (11) at arbitrarily 



A£ > IAS 



(11) 



large T. Thus in summary the MSF for chaotic dynamics 
has the same structure as for FPs and POs (the rotation 
symmetry follows from the rotation symmetry of each 
Ar(j)). ' 

Let us now discuss what the structure of the MSF 
means for the synchronizability of networks. We can clas- 
sify networks into three types depending on the magni- 
tude of the largest transversal eigenvalues 7 max in relation 
to the magnitude of the row sum a: (a) |7 max | < W\, (b) 
l7max| = M, and (c) |7 max | > \a\. 

As we have shown above, for stable synchronization 
it is necessary that |7 max | < r . If \j max \ > r , the syn- 
chronization is not stable. Since a is the eigenvalue of the 
coupling matrix associated with the synchronous mode, 
the MSF A max (a) describes the local stability within the 
SM, i.e., A max (cr) > for chaotic dynamics in the SM and 
Amax(f ) < for FPs or POs. This implies that |er| > ro 
in the first case and \a\ < r in the latter case. In other 
words, the row sum a gives an estimate of the critical 
radius ro. In particular, it allows us to give a complete 
classification (Table [I]). In networks of type (a) and (b) 
synchronization on a FP or a PO (stable within the SM) 
is always stable. For type (c) this dynamics may be sta- 
ble or not depending on the particular network topology 
(value of |7max|) and the dynamics in the SM (value of 
ro). On the other hand chaos synchronization is always 
unstable in networks of type (b) and (c) and it may be 
stable or not in networks of type (a) again depending on 
the particular network and the dynamics. 

Note that, in contrast to maps, autonomous flows with 
a stable PO in the SM always have ro = |oj, due to the 
PO's Goldstone mode. Thus for this case synchroniza- 
tion will be unstable for type (c) networks. For type (b) 



networks the stability of the synchronized solution in this 
case undergoes a destabilizing bifurcation. 

We now list some examples for the three types of net- 
works. The classification follows from the eigenvalue 
structure of the corresponding coupling matrices g^. 
Mean field coupled systems are of type (a), networks 
with only inhibitory or only excitatory connections are 
(up to the row sum factor) stochastic matrices and are 
thus of type (a) or (b) . Rings of uni-directionally coupled 
elements and two bidirectionally coupled elements are of 
type (b) and any network with zero row sum (a = 0) is of 
type (b) (trivial case) or (c) and therefore these systems 
can never exhibit chaos synchronization. 

Another conclusion we can draw from the structure of 
the MSF confirms the conjecture stated in [TD| : networks 
with er = are of type (c) and thus chaos synchronization 
is always unstable. 

Concerning the impact of noise on the delay-coupled 
network [16], for the case of FPs and POs stable syn- 
chronization will be robust to small noise strength. On 
the other hand, for the chaotic case there may exist an- 
other radius rb < ro, where the first UPO in the attractor 
loses its transverse stability and the attractor undergoes 
a bubbling bifurcation |17) fT8] . Then any network with 
r b < |7max| < r o will exhibit bubbling in the presence of 
small noise (or parameter mismatch), while any network 
with |7 ma x| < T\y will show stable synchronization, even 
in the presence of small noise. For large noise strength 
the linear theory cannot make predictions. 

Example - As an example we consider a network of 
optically coupled semiconductor lasers modeled by di- 
mensionless equations of Lang-Kobayashi [19 type 

E\t) = \{l + ia)n l {t) E l (t) + J2" =1 9lj&(t-T), 

Tn l (t)=p-n l (t)-(l + n l (t))\E l (t)\ 2 , (12) 

where E l and n l are the complex electric field ampli- 
tude and the inversion of the l-th laser, respectively. For 
our example, we choose the parameters as follows: Ratio 
between carrier and photon lifetime T = 200, injection 
current p — 10, a- factor a = 4. This results in a re- 
laxation oscillation period Tro ~ 28. Figure [TJd shows 
the A max = contour line of the corresponding MSF for 
networks with a — 0.4 for different values of the delay 
time r. For r = 20 (order of Tro) the contour line starts 
to become circular. For r > STro the shape of the MSF 
perfectly resembles our predictions. In this case we find 
ro < a = 0.4, i.e., the dynamics is chaotic. For r — 1 
and 10 the limit of large delay is not satisfied, hence the 
MSF does not exhibit the rotation symmetry. Note that 
for these values of the delay time the dynamics is a PO 
and since the system is a flow, the stability boundary 
reaches its maximum real value at re 1 ^ = \a\. 

Conclusion - We have shown that the MSF has a sim- 
ple universal structure in the limit of large delay: it is 
rotationally symmetric around the origin and either pos- 
itive and constant (if it is positive at the origin), or mono- 
tonically increasing and becoming positive at a critical ra- 



dius ro . This structure allows us to confirm a recent con- 
jecture |TD] about synchronizability of chaotic elements. 
Furthermore, we classify networks into three types de- 
pending on the magnitude of the maximum transversal 
eigenvalue of the coupling matrix in relation to the mag- 
nitude of the row sum. Importantly, this classification 
allows us to predict the synchronizability of general net- 
works of identical units with strongly delayed connections 
based solely on the modulus of the eigenvalues and the 
type of synchronized dynamics. In many cases this pre- 
diction is possible even without computing the critical 
radius ro (as shown in Table I). Although our results de- 
scribe the properties of coupled systems in the limit of 
large delay, practically they are expected to hold when 
the delay is two or three times larger than the character- 
istic timescale of the underlying system without delay. 
This is confirmed by our example as well as the results 
of Refs. [IT] , 
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